?getbbR Coding Workshop: 10th Meeting
GIS & Geospatial Data Analysis (Fall 2025)
R Coding Workshop for GIS: 10th Meeting
Outline
- Open Street Map
- Project checklist
- Raster Data
Open Street Map
- use place name to search for bounding box
taipei_bb_results <- getbb(place_name = 'taipei', format_out='sf_polygon')
taipei_bb_results| place_id | licence | osm_type | osm_id | lat | lon | class | type | place_rank | importance | addresstype | name | display_name | geometry |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 213710844 | Data © OpenStreetMap contributors, ODbL 1.0. http://osm.org/copyright | relation | 1293250 | 25.0375198 | 121.5636796 | boundary | administrative | 7 | 0.7009468 | city | 臺北市 | 臺北市, 臺灣 | POLYGON ((121.4571 25.10798… |
| 212454050 | Data © OpenStreetMap contributors, ODbL 1.0. http://osm.org/copyright | relation | 1527220 | 25.0119970 | 121.4656619 | boundary | administrative | 8 | 0.6192987 | city | 新北市 | 新北市, 臺灣 | MULTIPOLYGON (((121.2826 25… |
tpc_bb <- taipei_bb_results |> filter(place_id == 213710844)ntpc_bb <- taipei_bb_results |> filter(place_id == 212454050)
ntpc_bb |> mapview()osmdata
- create a query
assign("has_internet_via_proxy", TRUE, environment(curl::has_internet))tpc_mosque_results <- tpc_bb |>
opq() |>
add_osm_feature(key='building', value = 'mosque') |>
osmdata_sf()
tpc_mosque_resultsObject of class 'osmdata' with:
$bbox : relation(id:1293250)
$overpass_call : The call submitted to the overpass API
$meta : metadata including timestamp and version numbers
$osm_points : 'sf' Simple Features Collection with 67 points
$osm_lines : NULL
$osm_polygons : 'sf' Simple Features Collection with 3 polygons
$osm_multilines : NULL
$osm_multipolygons : 'sf' Simple Features Collection with 1 multipolygons
- leaflet visualization
leaflet() |>
addProviderTiles('CartoDB.Positron') |>
addPolygons(
data = tpc_mosque_results$osm_polygons,
color = "purple"
)tpc_mosque_polygon_sf <- tpc_mosque_results$osm_polygonsosmdata
kentucky_bb <- getbb('kentucky, USA')
kentucky_bb min max
x -89.57151 -81.96454
y 36.49712 39.14780
- Bounding Box
ky_results <- kentucky_bb |>
opq() |>
add_osm_feature(key = 'boundary', value = 'administrative') |>
add_osm_feature(key = 'admin_level', value = '4') |>
add_osm_feature(key = 'name', value = 'Kentucky') |>
osmdata_sf()
ky_resultsObject of class 'osmdata' with:
$bbox : 36.497118,-89.571509,39.1477997,-81.9645413
$overpass_call : The call submitted to the overpass API
$meta : metadata including timestamp and version numbers
$osm_points : 'sf' Simple Features Collection with 10411 points
$osm_lines : 'sf' Simple Features Collection with 266 linestrings
$osm_polygons : 'sf' Simple Features Collection with 0 polygons
$osm_multilines : NULL
$osm_multipolygons : 'sf' Simple Features Collection with 1 multipolygons
ky_sf <- ky_results$osm_multipolygons
ky_vars <- c("osm_id", "name", "ISO3166-2", "admin_level", "border_type",
"boundary", "nickname", "official_name", "ref", "ref:USCG", "ref:USPS", "ref:fips",
"short_name", "source:short_name", "type",
"wikidata", "wikipedia", "geometry")ky_sf <- ky_sf |> dplyr::select(all_of(ky_vars))
ky_sf |> mapview()ky_bookcase_results <- kentucky_bb |>
opq() |>
add_osm_feature(key='amenity', value = 'public_bookcase') |>
osmdata_sf()
ky_bookcase_resultsObject of class 'osmdata' with:
$bbox : 36.497118,-89.571509,39.1477997,-81.9645413
$overpass_call : The call submitted to the overpass API
$meta : metadata including timestamp and version numbers
$osm_points : 'sf' Simple Features Collection with 78 points
$osm_lines : NULL
$osm_polygons : 'sf' Simple Features Collection with 0 polygons
$osm_multilines : NULL
$osm_multipolygons : NULL
ky_bookcase_sf <- ky_bookcase_results$osm_points# ky_convex_sf <- ky_sf |> st_convex_hull() |> st_buffer()
# ky_convex_sf |> mapview()
ky_bbox_sf <- ky_sf |> st_bbox() |> st_as_sfc() |> st_as_sf(crs = 4326)
ky_bbox_sf |> mapview()Tessellation
library(tigris)To enable caching of data, set `options(tigris_use_cache = TRUE)`
in your R script or .Rprofile.
Attachement du package : 'tigris'
L'objet suivant est masqué depuis 'package:terra':
blocks
uac_sf <- urban_areas(year = 2020, progress_bar = FALSE)uac_sf |> st_crs() |> print()Coordinate Reference System:
User input: NAD83
wkt:
GEOGCRS["NAD83",
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4269]]
ky_sf |> st_crs() |> print()Coordinate Reference System:
User input: EPSG:4326
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
uac_sf <- uac_sf |> st_transform(4326)ky_uac_sf <- uac_sf[ky_sf,]
ky_urban_bookcase_sf <- ky_bookcase_sf[ky_uac_sf,]
ky_urban_bookcase_sf |> dim()[1] 64 26
mapview(ky_sf, col.region = 'lightblue', alpha.regions = 0.6, map.types = 'CartoDB.Positron') +
mapview(ky_uac_sf, col.regions = 'yellow') +
mapview(ky_urban_bookcase_sf, cex = 4)Observation Window
kybb_uac_sf <- uac_sf[ky_bbox_sf,]
kybb_uac_sf |> mapview()